Hydrogen molecule ion: Path integral Monte Carlo approach 
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Path integral Monte Carlo approach is used to study the coupled quantum dynamics of the electron 
and nuclei in hydrogen molecule ion. The coupling effects are demonstrated by comparing differences 
in adiabatic Born-Oppenheimer and non-adiabatic simulations, and inspecting projections of the 
full three-body dynamics onto adiabatic Born-Oppenheimer approximation. 

Coupling of electron and nuclear quantum dynamics is clearly seen. Nuclear pair correlation 
function is found to broaden by 0.040 ao and average bond length is larger by 0.056 oq. Also, 
non-adiabatic correction to the binding energy is found. Electronic distribution is affected less, and 
therefore, we could say that the adiabatic approximation is better for the electron than for the 
nuclei. 
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INTRODUCTION 

There is a number of phenomena in molecular and 
chemical physics which are influenced by the quantum 
behavior of both nuclei and electrons, rovibrational dy- 
namics being a good example, see Refs. [ll, 0, I3| and 
references therein. In case of light-mass nuclei, protons 
in particular, treatment of the quantum nature of the nu- 
clei is essential 0, H, |B| • This has proven to be important 
in the description of hydrogen bond, for example . 

Hydrogen molecule ion (Hj), being the simplest 
molecule, has been studied extensively [1| and it has often 
been used as an example or a test case for an improved 
method or accuracy In addition to the 

free molecule, influenced by an electric or rnagnetic 
field is a well-studied subject [li, [H [H [B, [H, [ll E^] . 
Furthermore, there is interest in descriptions that do not 
restrict to Born -Oppen heimer (BO) or other adiabatic 



approximations 
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Such extensions 



122,1231,12^ 

can be easily realized by using quantum Monte Carlo 
(QMC) methods [l^,!!!!, for example. 

Among the QMC methods the path integral formal- 
ism (PIMC) offers a finite-temperature approach together 
with a transparent tool to trace the correlations be- 
tween the particles involved. Though computationally 
extremely demanding, with some approximations it is ca- 
pable of treating low-dimensional systems, such as small 
molecules or clusters accurately enough. Some examples 
found in literature are H i2a], HD+ and H+ (s^, H2 clus- 

with special attention laid on 
140| . The approximations in these ap- 
proaches relate to the ad hoc type potentials describing 
the interactions between particles. 

In this work we evaluate the density matrix of the full 
three-body quantum dynamics in a stationary state and 
finite-temperature. This is what we call "all-quantum" 
(AQ) simulation. Secondly, the electronic part only 
is evaluated as a function of internuclear distance in 
the spirit of BO approximation, and thirdly, the adia- 
batic nuclear dynamics is evaluated in the BO potential 
curve. These allow us to demonstrate the non-adiabatic 




electron-nuclei coupling by a projection of the AQ dy- 
namics onto the adiabatic approximations. 

We need to approximate the — 1/r Coulomb potential 
of electron-nucleus interaction at short range to make 
calculations feasible. We realize this with a carefully 
tested pseudopotential (PP). Also, the absent (ortho) 
or negligible (para) exchange interaction of nuclei is not 
taken into account. Finally, we have to simulate a finite 
temperature mixed state. For convenience, we have cho- 
sen 300 K, but this essentially restricts the system to its 
electronic ground state. 

We begin with a brief introduction to the theory and 
methods in the next section. This includes description 
of the PP, and tools and concepts for the analysis in the 
following section. Then we carry on to the results. 



THEORY AND METHODS 

For a quantum many-body system in thermal equilib- 
rium the partition function contains all the information 
of the system 



41[. The local thermodynamical proper- 



ties, however, are included in the density matrix from 
which all the properties of the quantum system may be 
derived 4^ . The non-adiabatic effects are directly taken 
into account in PIMC. In addition, finite temperature 
and correlation effects are exactly included. 



Path integral Monte Carlo approach 

According to the Feynman formulation of the statis- 
tical quantum mechanics [i^ the partition function for 
interacting distinguishable particles is given by the trace 
of the density matrix. 
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where p{(i) — e~^^ , S is the action, /? = l/fceT, 
T = /9/M and Rm = Ro- M is called the Trotter num- 
ber and it characterizes the accuracy of the discretized 
path. In the limit M ^ oo we are ensured to get the 
correct partition function Z, but in practice sufficient 
convergence at some finite M is found, depending on the 
steepness of the Hamiltonian H. 

In the primitive approximation scheme of the PIMC 
formalism the action is written as 1441 
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+ U{R^,R,+ l■T), 
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where r) = ^[V{Ri) + V{R^+i)] and A = 

Sampling of the configuration space is carried out using 



the Metropolis procedure [45| with the bisection moves 
[4^. This way the kinetic part of the action is sampled 
accurately and only the interaction part is needed in the 
Metropolis algorithm. Level of the bisection sampling 
ranges from 3 to 6 in our simulations, respectively with 
the increase in the Trotter number. The bisection sam- 
pling turns out to be essential with large Trotter numbers 
to achieve feasible convergence, for nuclei in particular. 
Total energy is calculated using the virial estimator 47 1 . 
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Figure 1: (Color online) Hydrogen atom total energies with 
different Trotter numbers: infinite nuclear mass (triangle) 
and AQ (circle). Extrapolated ground state energies are 
-0.4947(f) Ha and -0.4938(3) Ha for infinite nuclear mass 
and AQ simulations, respectively. 



which in this work is replaced by a PP of the form [5C 
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Extrapolation of expectation values 

The Trotter scaling procedure [s^l for expectation val- 
ues is used to obtain estimates for energetics in the limit 
A/ — > 00. To use this procedure one needs expectation 
values with several different Trotter numbers. For the 
Trotter number Ad the scaling scheme is 



N 



(^)oo = {A)m + 



EC2i 



(3) 



where coefficients C2i are constants for a given tempera- 
ture and N represents the order of extrapolation. In this 
paper N = 2 has been used for the energies of H^, and 
= 3 for hydrogen atom energies, see Figs. [1] and [51 



Pseudopotential of the electron 



The parametres ac = 3.8638, a = 7.8857, a = 1.6617 
and h — —18.2913 were fitted using direct numerical so- 
lution to give the exact ground state energy of hydrogen 
atom and the wave function accurately outside a cut-off 
radius of about 0.6 cq. Also, a number of lowest energy 
orbitals of the hydrogen atom are obtained accurately 
outside the same cut-off radius 5l|. Because the bond 



length of H2 is about 2 oq, it is expected that bonding of 
the hydrogen molecule ion becomes properly described. 

Hydrogen atom reference energies for different Trot- 
ter numbers are shown in Fig[Tl where triangles are ob- 
tained from infinite nuclear mass and circles are from 
AQ simulations. Extrapolated ground state values are 
-0.4947(1) Ha and -0.4938(3) Ha for infinite nuclear 
mass and AQ simulations, respectively, statistical stan- 
dard error of mean (SEM) given as uncertainty in paren- 
thesis. We can note that within 2SEM limits proportion 
of these energies 0.9982 reproduces that of Rydberg con- 
stants, Rb./Roc = 0.9995. 



For the hydrogen molecule ion the potential energy is 

1 1 



^(ri,r2,R) 
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(4) 



where = |r — R,;|, i? = |Ri — R2I, r being the coor- 
dinates of the electron and R the internuclear distance. 
Eq. ^ sets challenges for PIMC arising from the sin- 
gularity of the attractive Coulomb interaction [i^ [4^, 



Spectroscopic constants 

Within the BO approximation of diatomic molecules 
the corrections to electronic energies due to rovibrational 
motion of the nuclei can be evaluated from a Dunham 
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polynomial [52| 

EvJ = - De + ^c{v +\)^ <^cXc[v + 

+ B,J{J +l)-a,J{J +l){v + ]^) + (6) 

where v and J are vibrational and rotational quantum 
numbers, respectively, and _Bc, ^cXc and are the 
spectroscopic constants. 

The spectroscopic constants of and are ob- 
tained as introduced in Ref. [s^ . In atomic units 
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Instead of determining these constants at the equilibrium 
distance only, as in Ref. [H^j , we evaluate expectation val- 
ues from the distribution of nuclei, e.g. for the rotational 
constant, 



(11) 



where the pair correlation function g{R) is normalized to 
unity. The other constants, Eqs. ([8|)- (fT0|l . are evaluated 
similarly. 



Centrifugal distortion 

Effects caused by the centrifugal distortion, arising 
from rotational motion of the nuclei, on the equilibrium 
distance can be assessed by inspecting the extremum val- 
ues of the energy of harmonic oscillator in rotational mo- 
tion: Ej{r) = ifc(r - rc)2 -l- J(J -l- l)/2fir'^. We find an 
approximate equation 
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4.B, 
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where i?o is the equilibrium distance. Eq. (jl2p . however, 
does not include the anharmonic effects shown in Eq. ([6|). 
which evidently increase the bond length. 

At finite temperature the rotational energy states 
should be weighted by the Boltzmann factor, which leads 
to 



AR = 



4^0 Ej-/(-/+l)exp(-/3BcJ(J+l)) 
M^c'^e' Ejexp(-/3SeJ(J-Kl)) 
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Figure 2: (Color online) potential curves with diflFerent 
Trotter numbers: M = 2" (square), M = 2^^ (triangle), 
M = 2^^ (circle), extrapolated values (dot), finite difference 
calculation with the pseudopotential (dashed) and with exact 
e~-p'^ potential (solid). 



where J = 0,1, 2, . . .. Using the spectroscopic constants 
from Ref. [54] , see Table HI and temperature of 300 K 
we obtain AR = 0.0043 ag. This approximation will be 
compared to our direct evaluation, below. 



RESULTS 

We consider three different cases separately in order 
to demonstrate the non-adiabatic effects. First, the elec- 
tronic part only is evaluated as a function of internuclear 
distance in the spirit of BO approximation. Secondly, the 
adiabatic nuclear dynamics is evaluated in the BO poten- 
tial curve. Finally, is treated fully non-adiabatically 
with the AQ simulation. These allow us to demonstrate 
the non-adiabatic electron-nuclei coupling by a projec- 
tion of the AQ dynamics onto the adiabatic approxima- 
tions. In addition, spectroscopic constants and isotope 
effects are looked into. 



Adiabatic electron dynamics 

Though the PP, Eq. (O , reproduces the hydrogen atom 
energy exactly, an error of —0.00342 Ha from the exact 
value —0.10263 Ha results in binding of another proton to 
form H J . This is demonstrated in Fig. [21 where potential 
curves of Hj from finite difference calculations with Vpp 
from Eq. ^ and exact V{r) = — are shown. 

Our PIMC energies with increasing Trotter number M 
and the extrapolation to M = 00 using Eq. ^ are shown 
in the same figure. These indicate clearly that the Trot- 
ter number has to be at least 2^^ in order to find the 
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Table I: Expectation values of spectroscopic constants, 
Eqs. (EJ-dlTJ. A Morse potential [s^l fitted to the FDpp 
potential curve is used in the evaluation of the energy deriva- 
tives. Corresponding pair correlation functions are shown in 
Fig. [3] First two columns are adiabatic nuclear dynamics re- 
sults and AQ results are in the third column. 







H+ (AQ) 


Ha cm^^ 


cm^^ 


cm^^ 



Be 



UJeXe 



0.0001366 
0.0001344 
0.0104816 
0.0104201 
0.0003552 
0.0003029 
6.445 X 10"' 
7.201 X 10"' 



30.35 15.24 
29.85705 

2328.96 1668.25 

2315.3 



78.92 
67.3 
1.432 
1.600 



35.33 



0.45 



29.26 

2229.77 
(2232)" 
90.73 

1.636 



This work 

Ref. [52] 
This work 

[52], [26]" 
This work 

Ref. [52] 
This work 

Ref. [52] 



minimum of the potential curve at the nuclear separa- 
tion R = 2.0 flQ. The extrapolated values are in good 
agreement with the potential curve FDpp, and there is 
almost a perfect match at i? = 2.0 oq, where the value of 
the extrapolated dissociation energy is 0.1061(2) Ha. 

For larger nuclear separations than 3.5 ao, however, we 
are not able to reproduce the potential curve with these 
Trotter numbers: we get too weakly binding molecule. 
This is assumed to be a consequence of the electronic 
wave function becoming more delocalized as the inter- 
nuclear distance increases, and thus the "polymer ring" 
representing the electron is not capable of sufficient sam- 
pling of configuration space. This error should diminish 
with increasing M. 

The electron-nucleus pair correlation function is shown 
in Fig. m and will be discussed below. 



Adiabatic nuclear dynamics 

For the quantum dynamics of the nuclei only (QN) we 
consider both and to see the isotope effect, too. 
The FDpp potential curve in Fig. [5] is used, for which 
convergence with respect to Trotter number is found at 
M > 2^ for both isotopes. Resulting pair correlation 
functions are shown in Fig. [3l 

Average nuclear separation of 2.019(1) aq for and 
2.007(2) ao for the isotope D^ is found with M > 2^. 
The full width at half maximum (FWHM) of the pair 
correlation functions are 0.539(1) ao and 0.454(1) ao for 
these isotopes, respectively. 

Difference in the bond length of between the adia- 
batic electron and adiabatic nuclei simulations, i.e total 
distortion, is 0.019 ao. Centrifugal contribution to this, 
the difference between one and three dimensional simu- 
lations of the nuclei, is 0.009(1) ao, which unexpectedly 
is about twice as much as the value 0.0043 ao evalu- 
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Figure 3: (Color online) Nuclear pair correlation functions: 
AQ (solid), Rt QN (dashed) and QN (dash-dotted). 
The difference in the average nuclear separation between QN 
and AQ is 0.056(3) ag. 



ated from the approximate Eq. (|13p . The anharmonic 
contribution, i.e. difference between total and centrifu- 
gal distortions, is 0.010(1) ao. 



In Ref. [54[ it was shown 
that anharmonic effects in H2 molecule contribute about 
the same amount to total distortion as centrifugal force, 
which turns out to be the case here, too. 

Difference between the total energies of the previous 
simulations (3D vs. ID) is 0.0009383(2) Ha, which is 
close to k^T « 0.00095 Ha as expected due to the pres- 
ence of the two rotational degrees of freedom in 3D. 
Difference between the dissociation energies of adiabatic 
electron and nuclear simulations, i.e. the zero-point vi- 
brational energy, is 0.0064(2) Ha. 

A Morse potential 53] fitted to the FDpp potential 
curve is used in the evaluation of the spectroscopic con- 
stants, see Table m This is justified because the nuclear 
simulations and analytical Morse wave function [H^ cal- 
culations coincide. The spectroscopic constants of Hj 
are close to those given in Ref. [H^] , which have been de- 
termined at the equilibrium distance of the nuclei, only. 
Same procedure is used for the spectroscopic constants 
of the other isotope. In Table[T]the same constants evalu- 
ated using the AQ instead of BO nuclear pair correlation 
function are also shown. 



Non-adiabatic "all-quantum" dynamics 



For H^ the total energy of AQ simulation with the 
Trotter number M = 2^^ is -0.60159(3) Ha. The extrap- 
olation procedure yields total energy —0.59872(3) Ha, 
which is only 0.0016 Ha more binding than the value 
—0.5971 Ha from variational Monte Carlo (VMC) simu- 
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Table II: energetics (atomic units). First three are BO 
and the next three are non-adiabatic values. 



Method 


-Etot 




^0 


R 


HF" 


-0.6026 


0.1026 




2.000 


VMC' 


-0.6026 


0.1026 




2.000 


PIMC" 


-0.6061(2) 


0.1061(2) 


0.0997(1) 


2.0 


VMC" 


-0.5971 




0.0971 


2.064 


MCDFT'' 


-0.581 




0.081 


2.08 


PIMC" 


-0.59872(3) 




0.09872(3) 2.075(2) 



" 5^: Hartree-Fock 

52|: VMC, Born— Oppenheimer 
" VMC, non-adiabatic 
"^0: MCDFT, non-adiabatic (SAO) 
"^This work 



lation [28] . The zero-point energy obtained from simula- 
tions is De - = 0.0074 Ha, see Table HH It should be 
pointed out that the error due to the pseudopotential in 
the AQ total energy is only about half of that found for 
the BO total energies. 

Difference in dissociation energies of AQ and the 3D 
QN simulations is 0.00097 Ha, which is about fceT 
revealing additional electronic energy degrees of freedom 
in the first. AQ simulation for HJ gives for the average 
nuclear separation R = 2.075(2) ag, which is 0.056 oq 
larger than that in the QN simulation. The AQ FWHM 
of the nuclear pair correlation function is 0.5785(2) ag, 
which shows a spreading of 0.040 oq compared to the QN 
results, see Fig. [31 

In Fig. |4] BO and AQ electron-nucleus pair correla- 
tion functions are compared. AQ projection onto the 
BO bond length, R — 2.0 aq, and BO results coincide, 
which indicates that the adiabatic BO approach for the 
electron dynamics is sufficient. Thus, it seems that the 
electron-nuclei coupling effects are more clearly seen in 
the dynamics of the nuclei, see Fig. [3l As one might ex- 
pect, there is a noticeable difference between the AQ and 
the BO electron-nucleus pair correlation functions due 
to varying bond length, see Fig. U) 

The AQ average nuclear separation is close to the value 
2.064 qq obtained by a non-adiabatic VMC simulation 
[2^. The AQ pair correlation function of the nuclei, see 
Fig. [21 coincides with the 5*^0 (Scaled Atomic Orbital) 
one in Ref. [2^ computed within the Multicomponent 
Density Functional Theory (MCDFT) scheme, not shown 
here. 

All the spectroscopic constants in Table [I] are defined 
using the derivatives from a fitted Morse potential, i.e. 
BO potential energy surface. Thus, the "AQ spectro- 
scopic constants" should be interpreted mainly as the di- 
rection of change in the values, except for Be- The expec- 
tation values of the spectroscopic constants are obtained 
by weighting the equations by the nuclear pair correlation 




Figure 4: (Color online) electron-nucleus pair correlation 
functions: AQ (solid, second lowest curve), AQ projection to 
R « 2.0 ao (sohd) and BO at i? = 2.0 ao (dashed). The 
latter two almost coincide. Dashed vertical line indicates the 
size of the pseudopotential core, r = 0.6 ao. For comparison 
corresponding pair correlation functions for hydrogen atom 
(dotted hue) and (dotted) obtained by using the analytical 
ground state wave function of hydrogen atom are also shown. 




Figure 5: (Color online) potential curves: Morse potential 
fitted to FDpp (dashed) and the elfective Morse potential ob- 
tained from the projection of the AQ simulation (solid), see 
the text for details. Corresponding nuclear pair correlation 
functions are shown in Fig. |3l The shift in the bond length is 
0.036 ao. 



function from the corresponding simulation. 

A projection of the AQ simulation to a potential curve 
of the nuclei is constructed with the help of the known 
solutions to the Morse potential. Distribution from the 
Morse wave function is fitted to the pair correlation 
function of the AQ simulation. The three-body sys- 
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Figure 6: (Color online) xy-plane (z-projection) snapshot 
from AQ simulation with Trotter number 2^^ for all particles. 
"Polymer ring" describing the electron is in the background 
and those of the nuclei are placed on top. 

tern is then presented by an effective two-body poten- 
tial. The projected potential curve shows clear differ- 
ences in the dynamics of the nuclei between BO and AQ 
simulations, see Fig. [5] The minima of the potentials 
are set to zero: the difference in the dissociation ener- 
gies between BO and the AQ projection is about 0.036 
Ha and the shift in the equilibrium distance is 0.036 oq. 
The spectroscopic constants with the projected poten- 
tial curve are Be = 29.26 cm"\ Uc = 2047.94 cm"\ 
WeXe = 78.12 cm~^ and ao = 2.110 cm~^. All this indi- 
cates that an effective Morse potential is not capable of 
describing non-adiabatic effects correctly. 

Finally, it may be of interest to see a visualization of 
the "polymer rings" representing the quantum particles 
in the PIMC simulation. So, Fig. [6] presents the xy-plane 
(z-projection) snapshot from AQ simulation with Trot- 
ter number 2^^ for all three particles. "Polymer ring" 
describing the electron is in the background and those of 
the nuclei are placed on top. 



CONCLUSIONS 

The three-body quantum system, hydrogen molecule 
ion (HJ), is revisited, once again. Path integral Monte 
Carlo (PIMC) method is used for evaluation of the sta- 
tionary state quantum dynamics. PIMC offers a finite- 
temperature approach together with a transparent tool to 
describe the correlations between the particles involved. 
We aim at tracing the electron-nuclei coupling effects 
in the three-body all-quantum (AQ), i.e. non-adiabatic, 
molecule. This is carried out by comparing the differ- 
ences in adiabatic Born-Oppenheimer (BO) and AQ sim- 
ulations, and inspecting the projections from the AQ sim- 



ulation onto the BO description of the electron-only and 
nuclear-only subsystems. 

The approach turns out to be computationally de- 
manding, but with the chosen pseudopotential for the at- 
tractive Coulomb potential and extrapolation to infinite 
Trotter number the task becomes feasible. By choosing 
low enough temperature, 300 K, we are able to compare 
our data to those from zero-Kelvin quantum methods 
available in literature. Among others we have evaluated 
spectroscopic constants and molecular deformation, also 
considering the isotope effects. 

With our fully basis set free, trial wave function free 
and model free approach we are not able to compete in 
accuracy with the zero-Kelvin benchmark values. How- 
ever, due to the mixed state density matrix formalism 
of PIMC we are able to present the most transparent 
description of the particle-particle correlations. 

Total energies from our simulations are more binding 
in nature compared to the benchmark values, see Ta- 
ble ini This is an expected effect of the pseudopotential 
in use, see Fig. [2] and FDpp therein. Quantum dynamics 
of the system is well described and distinct features of 
coupling are observed for the nuclei: shift of 0.056 ao in 
the equilibrium bond length, increase of 0.040 oq in the 
width of the pair correlation function of the nuclei and 
non-adiabatic correction of about 0.00097 Ha to disso- 
ciation energy. Electronic distribution, however, is less 
influenced by the coupling, see Fig. [51 and therefore, we 
could say that the adiabatic approximation is better for 
the electron than for the nuclei. 

Projection of the non-adiabatic three-body system 
with the help of Morse wave functions onto two-body 
nuclei-only subsystem indicates that Morse potential is 
not capable of describing non-adiabatic effects correctly, 
see Fig. [5] 
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